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We use the expansion-normalized variables approach to study the dynamics of a non-tilted Bianchi 
Type I cosmological model with both a pure, homogeneous source-free magnetic field and viscous 
fluid with both bulk and shear viscous effects, employing the ideal magnetohydrodynamic approx- 
imation. The dynamical system is then studied in detail through a fixed-point analysis which 
determines the local sink and source behaviour in combination with standard dynamical system 
theory and numerical experiments that give a very good indication of the asymptotic behaviour of 
the system, ft is further shown that for certain values of the bulk and shear viscosity and equation 
of state parameters, the model isotropizes at late times. 



I. INTRODUCTION 



^ ■ The current standard model of cosmology based on the Friedmann-LeMaitre- Robertson- Walker metric assumes that 

O^' the universe to a very high degree is spatially homogeneous and isotropic, which implies the existence of a perfect 
5_j I fluid matter distribution. However, if one wishes to formulate a cosmological model of the early universe, including 
bJO" anisotropic terms in the energy-momentum tensor is necessary. 

'~^, Viscous models have become of general interest in early-universe cosmologies largely in two contexts. Gr0n and 

p^ ■ Hervik (Chapter 13, [Ij) discuss these in some detail. The first relates through the idea of inflation through bulk 

^ viscosity. In models where bulk viscous terms are permitted to dominate, they drive the universe into a de Sitter-like 

fS| state. Because of these processes, the models isotropize indirectly through the massive expansion. Shear viscosity is 

^^ found to play an important role in universe models with dissipative fluids. The dissipative processes that result from 

^^ . shear viscous terms are thought to be highly effective during the early stages of the universe. In particular, neutrino 

^^ ' viscosity is considered to be one of the most important factors in the isotropization of our universe. 

^^ , Magnetic fields have also been thought to play a major role in the early universe, as it is well-known that after 

^^ ■ inflation, the early universe was a good conductor, even though the number density of free electrons dropped dramat- 

^^ [ ically during recombination, its residual value was enough to maintain high conductivity in baryonic matter. As a 

. . ' result, cosmic magnetic fields have remained frozen into the expanding baryonic fluid during most of their evolution 

J> ' (Page 115, 1^]). One can then analyze the magnetic effects on the dynamics of the early universe through the ideal 

magnetohydrodynamics (hereafter, referred to as MHD), in which the magnetic field source is considered to be a 

, ; perfect conductor, such that the energy momentum tensor is that for an ordinary magnetic field. 

C^ , Hughston and Jacobs [3| showed the in the case of a pure magnetic field, only Bianchi Types I, II, VI(/i = —1) 

(which is the same as Type III), and VII {h = 0) admit field components, whereas Types IV, V, VI {h = —1), VII 

{h y^ 0), VIII, and IX admit no field components. These results led to a number of papers of Bianchi models with a 

perfect-fiuid magnetic field source. 

With respect to the dynamical systems approach applied to the latter, which is what we are interested in studying in 
this work, LeBlanc studied Bianchi Type II magnetic cosmologies in which he provided an analysis on the future and 
past asymptotic states of the resulting dynamical system Q . LeBlanc also studied the asymptotic states of magnetic 
perfect-fiuid Bianchi Type I cosmologies [5]. Using phase plane analysis techniques, Collins studied the behaviour of a 
class of perfect-fluid anisotropic cosmological models, and established a correspondence between magnetic models of 
Bianchi Type I and perfect fluid models of Bianchi Type II [6,]. In addition, LeBlanc, Kerr, and Wainwright studied 
the asymptotic states of magnetic Bianchi Type VI cosmologies. In their work, they showed that there is a finite 
probability that an arbitrarily selected model will be close to isotropy during some time interval in its evolution [7]. 
We should note that Barrow, Maartens, and Tsagas did significant work in the reformulation of a 1-1-3 covariant 
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description of the magnetohydrodynamic equations that has provided further understanding and clarity on the role 
of large- vale electromagnetic fields in the perturbed Friedmann-LeMaitre- Robertson- Walker models [8;] . 

Viscous MHD Bianchi models, which are of interest to us have been presented in the literature on a number of 
occasions, van Leeuwen and Salvati studied the dynamics of general Bianchi class A models containing a magneto- 
viscous fluid and a large-scale magnetic field [9] . Banerjee and Sanyal presented some exact solutions of Bianchi Types 
I and III cosmological models consisting of a viscous fluid and axial magnetic field [iO| • Benton and Tupper studied 
Bianchi Type I models with a "powers-of-t" metric under the influence of a viscous fluid with a magnetic field [ll| . 
Salvati, Schelling, and van Leeuwen numerically analyzed the evolution of the Bianchi type I universe with a viscous 
fluid and large-scale magnetic field 12]. Ribeiro and Sanyal studied a Bianchi Type VIq viscous fluid cosmology 
with an axial magnetic field in which they obtained exact solutions to the Einstein field equations assuming linear 
relations among the square root of matter density and the shear and expansion scalars [13| . van Leeuwen, Miedema, 
and Wiersma proved that a nonrotating Bianchi model of class A containing a viscous fluid and magnetic field can 
only be of Type I or IVq [ij]. Pradhan and Pandey studied the Bianchi Type I model with a bulk viscous fiuid in 
addition to a varying cosmological constant [ll] . Pradhan and Singh studied the Bianchi Type I model in the presence 
of a magnetic field and shear and bulk viscosity, but assumed that the shear tensor was proportional to the expansion 
tensor [16] . Bali and Anjali studied a Bianchi Type I magnetized fiuid model with a bulk viscous string dust fiuid, in 
which they compared their results in the presence and absence of large-scale magnetic fields ^^ . 

We wish to study a Bianchi Type I non-tilted viscous magnetohydrodynamic model, but using the Hubble- 
normalizcd dynamical systems approach based upon the theory of orthonormal frames pioneered by Ellis and Mac- 
Callum [18|, which reduces the Einstein field equations, a coupled set of ten hyperbolic nonlinear partial differential 
equations to a system of autonomous nonlinear first-order ordinary differential equations. Such an approach has 
not been taken in the existing literature, which we believe will add valuable information to this already rich, but 
ever-expanding field. The advantage to using this approach is that one can obtain a very deep understanding of the 
evolution of these models at early times (near the initial singularity), at late times, and at intermediate times. This is 
mainly accomplished through a fixed-point analysis of the dynamical system where connections are then made with 
the global dynamics through sophisticated numerical experiments. Throughout, we assume that the signature of the 
metric tensor is (— , -I-, -I-, -I-), and the use of geometrized units, where G = c = 1. 

II. THE MATTER SOURCES 

For a viscous fluid magnetic cosmological model there are two contributions to the total energy-momentum tensor: 
the viscous contributions due to the shear and bulk viscosity and the magnetic field contribution. We will assume for 
simplicity that our fluid has negligible heat conduction, such that the heat flux vector vanishes. Then, the viscous 
fluid energy-momentum tensor is given by 

Vab = {^J■f + Pf)UaUb + gabPf " S^Hhab ~ ^rjaab, (1) 

where ^f,pf, and aab denote the energy density, pressure, and shear tensor due to the fluid, £, and rj denote the 
bulk and shear viscosity coefficients of the fluid, H denotes the Hubble parameter, and hab = UaUb + gab denotes 
the projection tensor corresponding to the aforementioned metric signature. A detailed derivation of Eq. ([T} can be 
found in (19.] . 

Following the arguments given on Pages 6 and 14 in [20|, we note that the energy- momentum tensor for an 
electromagnetic field is given by 

Tab = \uaUb{E^ + B^) + 2u^anlfucEgBd - EaEb - BaBb + hab (^' + B^) , (2) 

where n"-^'^'^ is the standard skew pseudo tensor, and Ea and Ba are the electric and magnetic field three-vectors 
respectively. Note that in an orthonormal frame where gab = nab == diag{—l, 1,1,1), the E^ and B^ terms in Eq. ©, 
take the form E^ = E'^Ea = E\ + E'^ + E'^, and B'^ = B'^Ba = -B^ + i?| -I- Bf. In addition, since we are assuming that 
the cosmological model is non-tilted, in both Eqs. ^ and ^, Ua is the four- velocity of a comoving observer such 
that u" = (1,0,0,0). 

For our purposes, we are to consider the ideal MHD approximation, where we assume that the early universe behaves 
as a perfect conductor. Then, the electric field, which is inversely proportional to the conductivity approaches zero, 
even if there is a current flowing. In other words, we are assuming that after recombination, the universe is such a 
good conductor that the cosmic electric fields required to drive current in it is negligible. Taking these arguments 
into account, Eq. ([2|) takes the form 

TBab = ^UaUb{B^)-BaBb + ^habB\ (3) 



The total energy-momentum tensor, denoted Tab, for our cosmofogical model is then 

Tab = Vab + Tsab, (4) 

where Vab is given by Eq. ^, and Tsab is given by Eq. ([3]). 

In order to formulation the evolution equations for our model, we will be required to know the total energy density, 
total pressure, and total anisotropic stress from Eq. ([4]), which we will respectively denote by /i,p, and iTab- Using 
the definitions 

/i = Ta6uV, p^h'^^'Tab, nab^KhtJcd-phab, (5) 



we find that 



A = M/ + ^ (Bf + Bl + Bl) , (6) 



p = wfif- 3^H +liB'i+ Bl + Bl) , (7) 



and 



iTab = -27]aab - BaBb + -hab {BJ + S| + B|) . (8) 



Note that in Eq. ([7]), we made use of the assumption that the fluid obeys the barotropic equation of state, pf — w^f, 
where — 1 < w < 1. 

Since we are dealing with expansion-normalized variables, we will re- write Eqs. ([S]) - (|H1) using the definitions [2l|: 

"-sl^' ^-WP- ft"-|J^ <^) 

We will also define the expansion-normalized magnetic field vector as 

B. = ^. (10) 

According to the definitions in Eq. (|9|), the expansion- normalized source variables take the form: 

n = nf + '^{Bl + Bi + Bl), (11) 

P^wrtf- 3^0 + ^ (S? + ^2 + Bl) , (12) 

and 

flab = -2m^ab - 9BaBb + 35ab {Bf + B^ + BI) . (13) 

In Eqs. ()12p and (|13p . ^o a-nd 770 are the expansion-normalized bulk and shear viscosity coefficients, and are defined 
according to 

1 = 3^0, f^Sryo. (14) 

Throughout this paper, we will assume that both ^0 and 770 are non-negative constants. In addition, in Eq. (J13p . Eab 
is the expansion-normalized shear tensor defined according to 'Eab = o'ab/H. 



III. BIANCHI TYPE I UNIVERSE DYNAMICS 

With the required energy- momentum tensor (Eq. Q), and the expansion- normalized source variables (Eqs. (jlip - 
P3p ) in hand, we will now derive the Bianchi Type I dynamical equations. The general evolution equations for any 
Bianchi type have already been derived in |21| and [22l | , and we will simply make use of their results in this section. 

The general evolution equations in the expansion-normalized variables using our notation are: 
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(16) 
As in Eq. ©, we have made use of the following notation: 

{^,,,R\N'\A,) = 1 {a,,,n\n'\a,) , (^,P,Q„ny) = ^ 0:*,^, g,,^,^) . (17) 

In the expansion-normalized approach, T^ab denotes the kinematic shear tensor, and describes the anisotropy in the 
Hubble flow, Ai and N'''^ describe the spatial curvature, while J7' describes the relative orientation of the shear and 
spatial curvature eigenframes. 

The dynamical system (|15p evolves according to a dimensionless time variable, r such that 

where H is the Hubble parameter. The deceleration parameter q is very important in the expansion-normalized 
approach, and through the evolution equation for H 

H' = -{l + q)H, (19) 

Eq. (1.90) in [23], and using Eq. (IT7|) . one can show that q is defined as 

q = 2j:^ + ^(n + 3P) 

= 2J:^ + nfl^l + ^-^yi^, + l{Bl + Bl + Bl), (20) 

where 2I]2 = i (S^jjE''''). 

In the case of a magnetic field source, one must also include an evolution equation for the magnetic field, which is 
the orthonormal frame analog of the standard Maxwell's equation. According to Eq. (71) in 24|, Eq. (2.4) in |7|, 
and Eqs. pU)) . p7|) . p^ . and p^ above, the magnetic field evolution is given by 

B', = Ba (-1 + q) + ^abB' + eab,Ii"B'' (21) 

The Bianchi Type I model is a flat anisotropic model and is Abelian, and therefore has the property that 

A' = 0, TVu - N22 - N33 - 0. (22) 

For convenience, we introduce the notation 

S+ = i (S22 + S33) , S_ - ^ (E22 - S33) , (23) 



such that E^ = E^ + E^. In the evolution equations ((T5|) . the expansion- normahzed angular velocity variables Ra 
can be found from the non-diagonal shear equations, ^,[2, 5^23 1 9'iid E'j^g. From these equations, we get that 

'^'--^^—' "^^-VSE^-SE/ "^^-VSE.+SE/ ^''^ 

To avoid situations where -\/3E_ — 3E+ = or a/3E_ + 3E+ = 0, we will set Bi = B3 = 0, and keep B2 ^ 0. Then, 
Ri = R2 — R3 — 0, and according to Eqs. (pi]) . (j^ . (f^ . and (PP)) . the evolution equations ([15]) and the constraint 
equations ([TC)) become: 

E; = -^62-fE+[g-2(l + ryo)], (25) 

E'_ = -^6| + E_[g-2(l + ,7o)], (26) 

6^ = 62 (-1 + 9 + \/3E_ + E+) , (27) 



where 



nj = 1- ^Bj - e2_ - E^ > 0, (28) 



q = 2{j:l + j:l) + nf(^ + ^\-^^o + lBl -i<w<i, eo>o, 7,0 > o. (29) 

The auxiliary equation, J7' in (J15p . after some algebra, takes the form 

n'f = nf{2q-l- 3w) + 4770 (E^ + E?.) + 9^0. (30) 

IV. A FIXED POINT ANALYSIS 

In this section, we consider the local stability of the equilibrium points of the system (P5|) - (l?f)l . The first thing to 
note is that the dynamical system Eqs. p5|)- (|27)) can be written as 

x' = f(x), (31) 

where x ~ [E+, E_, i52] G R-^, and f (x) denotes the right-hand-side of the dynamical system. 

The state space of the system is the subset of R^ defined by the inequality in Eq. (pS)) . which is equivalent to 

E^ + e2_ + ^Bl < 1, (32) 

so the state space is clearly bounded. This inequality also is a constraint for the initial conditions of the dynamical 
system. 

There is only one symmetry of the dynamical system and it is given by 

[E+,E_,62]^[S+,E_,-62], (33) 

so the system is invariant with respect to spatial inversions in the function B2, and we can therefore take B2 > 0- 
The critical points are points x = a that simultaneously satisfy 

f (a) = 0, (34) 

where f denotes the right-hand-side of the system (P5)) - (P7)) . The local stability is determined by linearizing this system 
at x = a, which leads to the relationship x' — Z?f (a)x. The stability of the system is then determined by finding 
the eigenvalues and eigenvectors of the derivative matrix Z?f(a). The restrictions on the parameters w,^o, and 770 are 
given by physically realistic limitations. Namely, we take the expansion- normalized viscosity coefficients ^Oi^o to be 
nonnegative constants. Any combinations of these parameters must additionally satisfy 11/ > 0, E_|_ e M, E_ G M, and 
/B2 > e M. In the work that follows, we will denote eigenvalues of the dynamical system by Ai, where i = 1, 2, 3, .... 



A. Equilibrium Point 1 



^ ^ l + 2rjo^ wjl + 2r]o) + yM-l + w) (9 + 4r/o - 12ry§ + 3w (-3 + 4770 + 47?§) + 36^o) 

+ 4(-l + w) 



^ ^ -3 - 67?o + ^;(3 + 67?o) + V^(-l + w) (9 + 47?o - 12r;g + 3w (-3 + 4r;o + 47yg) + 36go) 

4V3(-l + w) 
B2 = 0. (35) 

The cosmological parameters at this point take the form 

„ 4r?o + 9^0 on^ \ v2 3 - 3w + 47^0 + 9$o 

"^ = 3(-i+H ' '^ = 2(^+^°)' ^= — 3^3^;^ — ' ^^^^ 

where 

~l<w<l, 770 = 0, Co = 0. (37) 

With the restrictions ([37| . the equiUbrium point ((35|) takes the form: 



l + 3J(-l + w)^-w Vs (-1 + ^i^TTW + w) ^ ^ 

S- = ^ JT-, ^ -' ^2=0. (38) 



4(-l + w) ' 4(-l + w) 

In addition, the cosmological parameters p6p take the form: 

% = 0, q = 2, I]2 = 1, -l<w<l, (39) 

therefore, this point represents that of a Kasner circle, and we will denote it by /C. The eigenvalues of the dynamical 
system at this point are 

3 (-1 + ^/{^TTW + w) 

Ai=0, A2 = 3-3w;, A3 = ^ of , ^ ^ '■ (40) 

2(-l + w) 

We therefore see that in the region —l<w<l,£,o = riQ — 0,JC represents a local source, and can never be a local 
sink. We will subsequently denote the region defined by ([57)1 as U(K). 

B. Equilibrium Point 2 



E+ = 0, E_ = 0, B2 = 0. (41) 

The cosmological parameters at this point take the form 

% = 1, q^^{l + 'Sw~9^o), S2 = 0. (42) 

Clearly, this point represents the flat FLRW universe, which we will denote by JF. The eigenvalues of the dynamical 
system at this point are given by 

Ai = i(-l + 3u;-9eo), A2 = A3 = -(-3 + 3u;-477o-9eo), (43) 

where in Eqs. (gU and gS]), 770 > 0, ^o > 0, and -1 < u; < 1. 
The point J^ represents a local sink if 

% > 0, Co > 0, -1 < u; < -, (44) 



%>0, ^<^<1, ^o>^(-l + 3u;). (45) 

In Fig. ([!]), we have denoted the region defined by (l44l) and (|45)) as S1{F). 

An important point to note is that q = —1 when < ^o < § and w = 3^o ^ 1, thus, the equihbrium point in the 
domain defined by these values of rjo, ^o jand w does not correspond to a self-similar solution. In particular, if one 
chooses ^0 = such that w = —1, the corresponding model is locally the de Sitter solution [2]. 

The point J" represents a saddle point if 

m = 0, ^<w<l, 0<eo<^(-l + 3w;), (46) 

or 

7/0 = 0, w = l, 0<^o<^, (47) 

or 

%>0, ^<w<l, 0<eo<^(-l + 3u;), (48) 

where in each case Ai > and A2 = A3 < 0. We will subsequently denote the region defined by (|15)) - (H5)) as SA(F). 
The point J^ can also represent a local source if 

Vn^O, ^" = 1, Co = 0, (49) 

where in this case, Ai > and A2 = A3 = 0. IWe will subsequently denote the region defined by Eq. (|49l) as U(F). 

C. Equilibrium Point 3 



We will denote this equilibrium point as BI_ 



MV- 



E+ = - ! fo^a ^^ f 3^ ~ 52770 + 12772 + 9{w + 2w7jof + w {36 + 48770 - 48772) 

16(5 - 6770 + ^(3 + 67/0)) V 

+ y/(5 - 6?7o + 777(3 + 6?7o))2 (81 - 67i;(3 - 2770)2 - 28770 + 477^ + 9{w + 2wr]o)^ + 288^0) ) , 



E_ = - , ,„ , ^ .. VS 35 - 52770 + Wo + 9{w + 2wm? +w{36 + 48770 - 487/2) 

16(5 - 6770 + w{3 + 6770)) V 

+ ^^(5 - 6770 + 77;(3 + 677o))2 (81 - 67x;(3 - 2770)2 - 287/0 + 47/2 + 9{w + 2wr]o)^ + 288^0) j , 

B2 = ^[V ( -51 + 127«(1 - 27/o)2 + 52770 - 127/2 _ g(^ _^ 2it;77o)2 - 144^0 

4V3 V V (50) 

- ^^(5 - 6770 + w;(3 + 677o))2 (81 - 677;(3 - 2r/o)2 - 28770 + 4772 + 9{w + 2wtjoY + 288^o) ) ) • 

The cosmological parameters at this point take the form 

"/ = - ,.., , 1 t^^a ^^ f-45 + 114% - 82772 + 24773 - 1277;(1 - 2r7o)2(l + 2770) + 977;2(l + 2770)^ 
16(5 - 6770 + w[3 + 6770)) V 

+ \/(5 - 6770 + 77;(3 + 677o))2 (81 - 6w(3 - 2770)2 - 28770 + 4772 + 9(w + 2wr]of + 288Co) 
+ 2770^/(5 - 6770 + 77;(3 + 6r/o))2 (81 - 67«(3 - 2770)2 - 28?7o + 477^ + 9(w + 27^770)2 + 288Co) j , 



and 



where 



1 



9 = 



4(5 - 67?o + w(3 + 6m)) 



55 - 76r;o + Urj^ + 9(w + 2wt]o)^ + w (48 + 72r;o - 4877^) 



+ \/(5 - 6770 + w(3 + 677o))2 (81 - 6w{3 ~ 2770)2 - 28r;o + A'qi + 9{w + 2wria)^ + 288^0) ) , 



1 



64(5 - 6770 + u;(3 + 677o))2 



35 - 52770 + 12r;^ + 9(w + 2u;77o)^ + w (36 + 48770 - 48r;^) 



3 . n 1 

'yo > 2' Co = 0, w = -, 



(51) 



+ y^(5 - 6r/o + w{3 + 6?7o))2 (81 - 6w;(3 - 2770)2 - 28% + 477^ + 9{w + 2wm)^ + 288^0) j , 



(52) 



3 1 -5 + 6770 

'^°>2' 3<"<TT6;^' 



0<eo< g(-l + 37i;) 



(53) 



We will subsequently denote the region defined by ^, ^ as S2(BIMV). 

Under the conditions (j52p . the characteristic polynomial of the dynamical system at this point is given by 



c(A) 



^-^ (-36 (3 + v/(-3 + 2r/o)2) - 877i^(-4 + A) + (87 + 28v/(-3 + 2r/o)2) A - 2 (9 + 4^(^3 + 2770)2) A^ 

+ 3A=' - 4r/2 (4 (9 + V(-3 + 2770)2) - 13A + 2A2) 



-3 + 2770 



- 2770 (-12 (9 + 2V(-3 + 2770)2 + (59 + 12V(-3 + 2770)2 A - 12A2 + A 



Solving Eq. ((54|) for A in the region described by Eq. ([52|). gives the eigenvalues 

Ai-0,A2--l-2r;o, 



(54) 



(55) 



which shows that in this case, the equilibrium point is a local sink. 

We were not able to find closed-form expressions for the eigenvalues of the dynamical system in the region defined 
by Eq. ([53)1 . We therefore used numerical methods to determine some eigenvalues to get a good idea as to the local 
behaviour of the system about this point. As can be seen from Table I, all of the computed eigenvalues in this region 



TABLE I: Eigenvalues of the Dynamical System determined by numerical methods in the region defined by Eq. (f53l 



Ai 


A2 


A3 


-6.19013 


-6.17073 


-0.0200152 


-8.19018 


-8.17053 


-0.0199997 


-10.1902 


-10.1707 


-0.199906 


-12.1903 


-12.1708 


-0.0199846 


-14.1903 


-14.171 


-0.0199804 


-16.1904 


-16.1711 


-0.0199773 


-18.1904 


-18.1713 


-0.0199751 


-20.1905 


-20.1714 


-0.0199733 



are negative. There is therefore strong evidence that the equilibrium point in this region is that of a local sink. 

It is of interest to note that to the authors' best knowledge at the time of the writing of this paper, the equilibrium 
point of our dynamical system with S+7^0, S_7^0, and 62 7^ has not been given previously in the literature in 
the context of a dynamical systems approach to modelling viscous magnetic fluid universes. This in turn also implies 
that the model will not isotropize in the regions defined by Eqs. (|52)) and (|53|l. 



V. BIFURCATIONS 

The physical equihbria found in the previous section are related to each other by a sequence of bifurcations, which 
can be understood as follows. Consider the linearized equation for B2 at J^ 



B', = B2 



i-l + Sw-9^0) 



(56) 



which destabilizes T at 



^o=g(-l + 3u;). 



(57) 

We therefore see that the line given by Eq. (1571) is very important with respect to the dynamical evolution of the 
system as it governs the bifurcations and transitions from one equilibrium state to the next. Interestingly, the quantity 
given by Eq. (|57|l is a boundary for the regions defined by each equilibrium point, J" and BXmv- In the case of /C, 
the bifurcation equation is trivially satisfied for ^0 = 0, and w = ^ which satisfies the boundary region for K.. 

For convenience, we have summarized the results of the previous section in Fig. ([1]) which depicts the different 
regions of sinks, saddles, and sources of the dynamical system. 



FIG. 1: A depiction of the different regions of sinks, sources, and saddles of the dynamical system as defined by the aforemen- 
tioned restrictions on the values of the expansion-normalized bulk and shear viscosities, ^o, ^70 and equation of state parameter, 
w. 
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As can be seen from Fig. ([IJ and the above fixed-point analysis, a possible bifurcation sequence is given for when 
i < w < 1 is held constant, and % > |. Then, we have that 

BIMV ^^^^^ ^ BIMV ^^^-^-^ 4 Sl{F). 

Therefore, the dynamical system changes its equilibrium behaviour as the expansion-normalized bulk viscosity coef- 
ficient, ^0 changes, which in turn depends on the chosen value of w, in addition to the requirement that 770 > |. 

VI. QUALITATIVE PROPERTIES OF THE SYSTEM 

An important question to ask in analyzing some qualitative properties of the dynamical system is whether there 
are any invariant sets of the dynamical system. A very useful proposition in this regard is given by Proposition 4.1 
in [23], which states that for a dynamical system of type ([5T|). if Z : R" ^ R is a C^ function such that Z' — aZ, 
where a : R" — > R is a continuous function, then the subsets of R" defined hy Z > 0, Z ~ 0, and Z < are invariant 
sets of the flow of the system of differential equations. Looking at Eq. (P7|) . we see that this proposition applies with 
Z = B2, so B2 — and B2 > arc invariant sets of the system. We also note that if one sets rjQ — (^0 — in Eq. (|30p . 
then the proposition also applies with Z — flf, with flf > 0, defining invariant sets of the system. 

With respect to existence of limit sets and monotone functions, we first make a proposition about the late-time 
dynamics of the system. 

Proposition 1 Consider the dynamical system US]) in the region defined by —1 < w < ^, £,0 — ^> o.nd rjQ > 0. Then, 
as T ^>- -fcx), S^ = S^ + S'L — > and y8| — > 0, hence, the model isotropizes. 

Proof. The details of the proof essentially follow the arguments given in the appendix of [23 . Looking at Eq. (|29|) , 
we have upon using Eq. (j28p that 

^2/3 — 3?i;\ 2^3 — 9w"\ , 3w + 1 9. 



^'{^ +^H^ +^--2^0. (58) 



The ft'f evolution equation (15(71) now takes the form 



n'j ^rif 



j:-^i3-Sw)+BU^^]~9£o 



4770S]" + 9^0. (59) 



2 
In addition, from the generalized Friedmann equation, Eq. (|28p . we have that 

nf = 1- \bI y? 

I 

^ fi/ - 1 < 

^VLf < \. (60) 

Therefore, to prove the proposition all we have to do is show that the function 51j is monotonically increasing, that 
is, XI f > 0. Since we have the freedom to choose ^0, we will let ^0 = 0. Then, 



^%-l = -^^2-^2 



S^ (3 - 3w) + Bl' 



-4?7oE2 > 0^ -1 < w < -. (61) 

o 

Therefore, in the region where 770 > 0, ^0 = 0, and — 1 < w < -I, 51^ is monotonically increasing. In can therefore be 
said that in this region, 

lim Q.f = 1. (62) 

r— >-+oo 

Using this fact then in the Friedmann equation, (pS)) . we have that 

lim % = 1^ lim ( — BI-tA=Q. (63) 
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The latter then impUes that 

hm E^ ^ lijn ^2 ^ Q ^Q^^ 

To be technically correct, one should really set 770 = 0. This is because we have just shown that in the region under 
consideration, the model isotropizes. Therefore, the isotropic universe necessarily has a perfect fluid distribution, 
hence 770 = 0, which is still perfectly valid under our assumptions. ■ 

Note that, the region under consideration of the proposition is precisely one of the regions for which S1{F) is a 
local sink, which is precisely the conclusion we reached in the proof of the proposition above. 

One can also use the extended LaSalle invariance principle, that is extending the LaSalle invariance principle for 
negatively invariant sets to get an idea of the past asymptotic behaviour of the dynamical system. The form of the 
LaSalle invariance principle that we will require for this purpose is given as Proposition B.3. in [26|. Namely, consider 
an autonomous system of first-order differential equations x' = f(x), and let Z : M" — ?> M be a C^ function. Let 
S C M" be a closed, bounded, and negatively invariant set, and Z'{x.) = VZ ■ f(x) < 0, V x e 5, then V a e 5, 
a(a) C {x e S\Z'{x) = 0}. 

Proposition 2 For the dynamical system 11 15]). a(a) C {flf = 0}. 

Proof. Clearly the set {ftf = 0} is negatively invariant, closed, and bounded. Eq. ([5^1 in the region defined by 
{ilf — 0} reduces to 4ryoS'^ + 9^o- Clearly, this function is < if and only if ryo = Co = 0. Therefore, fl'f = if 
^f — £,0 = Vo — 0; which is precisely the region defining the Kasner circle according to Eqs. pB]) and Eqs. (|39|) . and 
so a(a) C {flf = 0}. ■ 

Because the dynamical system (jlSp is a three-parameter family of autonomous ordinary differential equations, hence 
having complex behaviour, we were not able to provide much more of a qualitative analysis beyond what we have done 
so far with regards to asymptotic behaviour and the computations of the local sinks and sources. For the remainder 
of the paper, we therefore employ numerical techniques to study the dynamical system more thoroughly to obtain a 
solid understanding of the complex dynamics involved as related to the three parameters. 

VII. A NUMERICAL ANALYSIS 

The goal of this section is to complement the preceding stability analysis of the equilibrium points with extensive 
numerical experiments in order to confirm that the local results are in fact global in nature. For each numerical 
solution, we chose the initial conditions such that the constraint Eq. ([28|) in addition to ^2 > were satisfied. 

Although numerical integrations were done from < r < 3000, for demonstration purposes we presented solutions 
for shorter time intervals. We completed numerical integrations of the dynamical system for physically interesting 
cases of u; equal to (dust), 0.325 (a dust/radiation mixture), and 1/3 (radiation). Also note that in the subsequent 
plots, asterisks denote initial conditions. 
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A. Dust Models: w = 



1. Co = 0.1,r?o = 0.2 



FIG. 2: This figure shows the dynamical system behaviour for ^o = 0.1, 770 = 0.2, and w = Q. The diamond indicates the 
FLRW equihbrium point, and this numerical solution shows that it is a local sink of the dynamical system. The model also 
isotropizes as can be seen from the last figure, where E±, B2 — >■ as r ^^ 00. 
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2. ^o = l,??o = 0.5 



FIG. 3: This figure shows the dynamical system behaviour for ^o = 1, ?7o = 0.5, and ui = 0. The diamond indicates the FLRW 
equilibrium point, and this numerical solution shows that it is a local sink of the dynamical system. The model also isotropizes 
as can be seen from the last figure, where E±, S2 — ^ as r — !> cxd. 
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B. Radiation Models: w = 1/3 

1. Co = 1.5, 770 = 



FIG. 4: This figure shows the dynamical system behaviour for ^o = 1-5, 770 = 0, and w = 1/3. The diamond indicates the 
FLRW equihbrium point, and this numerical solution shows that it is a local sink of the dynamical system. The model also 
isotropizes as can be seen from the last figure, where E±, S2 — !■ as r ^^ 00. 
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2. Co = 1.5,r?o=0.5 



FIG. 5: This figure shows the dynamical system behaviour for ^o = 1-5, 770 = 0.5, and w = 1/3. The diamond indicates the 
FLRW equilibrium point, and this numerical solution shows that it is a local sink of the dynamical system. The model also 
isotropizes as can be seen from the last figure, where E±, B2 — > as r -> 00. 
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Co = 0,r?o = 2 



FIG. 6: This figure shows the dynamical system behaviour for ^o = 0, t^o = 2, and w — 1/3. The circle indicates the BIMV 
equilibrium point, and this numerical solution shows that it is a local sink of the dynamical system. The model does not 
isotropize with respect to the anisotropic magnetic field as can be seen from the last figure, where £2 > as r — >■ 00, but does 
isotropize with respect to the spatial anisotropic variables, Ej-,^- as r — >■ 00. This state is also special, since according to 
our fixed-point analysis, this behaviour is only exhibited for w — 1/3, -qo > 3/2, and ^0 = 




"r- 
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I Co = 0, r?o = 10 



FIG. 7; This figure shows the dynamical system behaviour for ^o ~ 0, J70 = 10, and w = 1/3. The circle indicates the BIMV 
equilibrium point, and this numerical solution shows that it is a local sink of the dynamical system. The model does not 
isotropize with respect to the anisotropic magnetic field as can be seen from the last figure, where £2 > as r — >■ 00, but does 
isotropize with respect to the spatial anisotropic variables, E±, ^- as r — >■ 00. This state is also special, since according to 
our fixed-point analysis, this behaviour is only exhibited for w = 1/3, -qo > 3/2, and ^0 = 
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C. Dust/Radiation Models: w = 0.325 
1. Co = 0.5, r?o = 0.5 



FIG. 8: This figure shows the dynamical system behaviour for ^o = 0.5, rjo = 0.5, and w = 0.325. The diamond indicates the 
FLRW equilibrium point, and this numerical solution shows that it is a local sink of the dynamical system. The model also 
isotropizes as can be seen from the last figure, where E±, B2 — i- as r ^^ 00. 
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D. Dust/Radiation Models: w = 0.325 



FIG. 9: This figure shows the dynamical system behaviour for ^o ~ l, Vo = 2, and w = 0.325. The diamond indicates the 
FLRW equilibrium point, and this numerical solution shows that it is a local sink of the dynamical system. The model also 
isotropizes as can be seen from the last figure, where E±, S2 — !■ as r — !> 00. 
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E. Bifurcations 



As we alluded to previously, according to Eq. ((57)) . the line Co = | (^1 + Sw) governs bifurcations of the dynamical 
system. We therefore thought it would be constructive to display this bifurcation behaviour for cases where first 
Co < 5(— l + Sic), then Co = ^(—1 + 3^;), and finally, Co > 5(— l + 3w). For the purposes of this numerical 
experiment, we specifically chose w = 1/2, 770 = 5, therefore requiring that the three aforementioned cases reduce to 
Co < igi Co = igi and Co > xg- 



FIG. 10: These figures show bifurcation behaviour for a varying expansion-normalized bulk viscosity coefficient, ^0, while w 
and ?7o were held fixed. The circles indicate BIMV equilibrium points, while the diamond indicates the FLRW equilibrium 
point. For the first figure, ^0 ~ 0.05, for the second figure, Co = jg, and for the last figure, Co = 0.6. Note how the increasing 
values of Co first result in a slight shift of the equilibrium point position of BIMV, and then finally a transition to a new state, 
namely the FLRW equilibrium, which was predicted by our fixed-point analysis. 
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VIII. CONCLUSIONS 

We have presented in this paper a comprehensive analysis of the dynamical behaviour of a Bianchi Type I viscous 
magnetohydrodynamic cosmology, using a variety of techniques ranging from a fixed point analysis to analyzing 
asymptotic behaviour using standard dynamical systems theory combined with numerical experiments. We were able 
to show that for cases in which ?7o > 0, Co > 0, —1<w<^ott]q>0, ^ < w < 1, Co > | (—1 + 3w), the 
dynamical model isotropizes asymptotically, that is, the spatial anisotropy and the anisotropic magnetic field decay to 
negligible values giving a close approximation to the present-day universe. We were also able to show that for regions 
in which 770 > |, Co = 0, w = ^ or 770 > |, ^ < w < ~^^''° , < Co < | (— 1 + 3w), the model does not isotropize, 
rather at late times goes into a stable equilibrium in which there is a non-zero magnetic field. 

Further to our analysis, the flat FLRW model which was our S1(F) equilibrium point, is of primary importance with 
respect to observational/physical significance. Through our fixed point analysis, we showed that S1(F) represents a 
saddle point if 770 = 0, ^ < w < 1, < Co < i(~l + 3w), tjq — 0, w = 1, < Co < |, or 770 > 0, ^ < 
u! < I, < Co < I (— lH-3w), (which was denoted above by SA(F)). In these regions, S1(F) attracts along its 
stable manifold and repels along its unstable manifold. Therefore, some orbits will have an initial attraction to S1(F), 
but will eventually be repelled by it. There is therefore a time period during which the cosmological model will 
asymptotically isotropize, and be compatible with present-day observations of high-degree isotropy. 
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